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ABSTRACT 


Prony* s  method  is  a  simple  procedure  for  determining  the  values  of  para¬ 
meters  of  a  linear  combination  of  exponential  functions.  Until  recently, 
even  the  modern  variants  of  this  method  have  performed  poorly  in  the  presence 
of  noise.  We  propose  a  simple  procedure  for  estimation  of  the  signal  para¬ 
meters  in  the  presence  of  noise.  This  procedure  is  very  close  in  form  and 
assumptions  to  Prony's  method.  However,  in  preliminary  tests,  the  performance 
of  the  method  ie  close  to  that  of  the  best  available,  more  complicated,  approaches 
which  are  based  on  maximus  likelihood  or  on  the  use  of  eigenvector  or  singular 


vector  decompositions. 


Z.  Introduction 

Nearly  two  hundred  yeara  ago  Prony  developed  a  simple  procedure  for  de¬ 
termining  the  values  of  parameters  of  a  linear  combination  of  exponential 
functions  fromunifoxmly  spaced  samples  £l].  Today  "Prony1 s  method"  is 
usually  taken  to  mean  the  least  squares  extension  of  the  method  as  presented , 
for  example,  by  Hildebrand  £2!. 

A  short  record  of  a  data  sequence  y(n),  n  ■  1,2,...N,  is  assumed  to 
be  composed  of  uniformly  spaced  samples  of  a  sum  of  exponential  signals,  x(n), 


and  measurement  noise  w(n).  That  is, 

y(n)  •  x(n)  +  w(n)  for  n  *  1,2,...,N  (1) 

where 

x(n)  -  Z  a(k)(C(k))“  w 

k-1 

N  >  2M  (3) 

C(k)  ■  exp(s(k))  (4) 

The  values  of  the  signal  parameters  a(k)  and  s(k)  for  K  ■  1,2 . M  are 


unknown  complex  numbers.  Often  the  value  of  M  is  unknown  also.  However,  let 
us  initially  assume  that  the  value  of  M  is  known. 

Following  the  derivation  of  Hildebrand  £2]  we  note  that  the  signal  x(n) 
satisfies  a  linear,  homogeneous  difference  equation  with  constant  coefficients 

M 

Z  b(k)  x(n-k)  ■  0  for  M  <  n  <  N  (5) 

k«0 


where 


b(0)  -  1  (6) 


The  roots  of  the  prediction-error-filter  polynomial  B(s)  provide  the  values  of 


(7) 


the  exponent  parameters  C(k),  and  hence  s(k): 

M  -k 

B(z)  -  Z  b(k)a 

k-0 

- ;  rw*».-n 

k-1 

Hildebrand  explicitly  considers  noisy  data  and  specifies  Prony's  method 
by  the  following  three  steps: 

(1)  Using  the  method  of  least  squares,  minimize  the  approximation  error 

N  M  A  2 

A  *  Z  I  b(k)  y(n-k)  +  y(n)  (8) 

n-M+1  k-1 

by  best  choice  of  the  coefficients  b(k)  {_  3^*  For  N  >  2M  and  for  noisy  data 
the  solution  will  be  unique  with  high  probability.  However,  if  the  resulting 
set  of  normal  equations  is  singular,  then  the  pseudo  inverse  of  the  coefficient 
matrix  can  be  used  to  choose  the  minimum  norm  solution. 

(2)  After  the  M  values  of  £(k)  are  determined,  the  roots  of  the  pre- 


diction-error-filter  polynomial,  B(Z),  are  found,  using  £(0) 


M  -  -k 
Z  b(k)zk 

k-0 


(1  -  C(k)*"1) 


The  corresponding  exponent  values,  s(k),  can  then  be  found  from  formula  (4). 

(3)  Having  determined  the  values  C(k)  for  k-  1,2, ...,M,  the  error  in 

approximating  the  observed  data  by  a  linear  combination  of  exponential  signal 

components  then  becomes  linear  in  the  M  values  of  a(k): 

M  n 

e(n)  -  y(n)  -  Z  «(k)  fc(k)l  for  n  -  1,2,... N  (10) 

k-1  L  J 

The  M  estimates,  a(k),  can  be  determined  by  minimizing  the  summed, 


magnitude-squared  error: 


4 


E  » 


N 

I 

n-1 


y(n)  - 


M 

l  a 

k-1 


(k)  [c(kjj! 


(11) 


It  it  veil  known  that  the  errors  in  signal  parameters  which  are  estimated 
by  Prony'a  method  can  be  discouragingly  large  £2, 4.]  •  For  insight  into  this 
phenomenon  we  recanaend  calculation  and  study  of  the  Cramer-Rao  (CR)  bounds 
for  the  variance  of  the  error  in  the  estimated  parameters  and  comparison  of 
the  threshold  of  estimation  of  the  Prony  method  with  that  of  the  maximum  like¬ 
lihood  method.  C.7,8,9,loJ.  By  the  threshold  of  estimation,  we  mean  the  value 
of  signal-to-noise  ratio  (SNR)  at  which  the  variance  of  an  estimation  error 
begins  to  depart  very  rapidly  from  the  corresponding  CR-bound  value. 

As  another  example  of  the  application  of  Prony'a  method  consider  the 
problem  of  estimating  the  parameters  of  a  zero-mean,  autoregressive,  moving 
average(ARMA)  stationary  random  sequence  from  estimates  of  its  covariance 
values.  Various  investigators  have  recognized  that,  after  a  finite  number  of 
lags,  the  true  underlying  covariance  values  satisfy  a  linear,  homogeneous, 
difference  equation  with  constant  coefficients  |  11,12,13].  That  is,  after 
a  finite  number  of  lags,  the  estimated  covariance  sequence  can  be  represented 
as  a  linear  combination  of  exponentials  (i.e.  the  true,  underlying  covariance 
sequence  which  satisfies  the  homogeneous  difference  equation)  plus  measure¬ 
ment  noise.  This  measurement  noise  may  be  only  the  error  sequence  in  estimating 
the  covariance  values  from  a  finite  observation  of  the  ARMA  sequence.  Part 
of  it  might  also  be  due  to  additive  noise  in  the  observation  of  the  ARMA 
sequence.  Some  of  these  ideas  have  been  restated  by  Cadzow  fl4,  15*] . 


II.  Prony  Methods  for  Noisv  Data 

In  previous  and  related  work,  we  have  shown  how  one  can  extend  the  threshold 
of  estimation  of  Prony' •  method  to  much  lower  values  of  SNR  and  how  one  can 


improve  parameter  estimation  at  values  of  SNR  above  this  threshold^, 9, 10,11, 16-22] 
The  major  source  for  these  improvements  is  the  use  of  information  about  the 
rank,  M,  of  a  matrix  of  signal  covariance  values  or  a  matrix  of  samples  of 
the  signal.  If  there  is  no  prior  information  about  this  rank,  it  is  estimated 
from  the  data  using  singular  value  decomposition  (SVD).  The  most  important 
computational  step  is  a  preprocessing  step,  before  application  of  Prony's 
method.  A  prediction  order  L  which  is  larger  than  the  value  of  M  is  chosen. 

The  measured  data  matrix  or  the  matrix  of  estimated  covariance  values  is  re¬ 
placed  by  a  matrix  of  the  prescribed  rank,  M,  which  is  the  best  least  squares 
approximation  to  the  given  matrix.  Other  investigators  have  presented  closely 
related  approaches  [23-29]. 

In  this  work  we  advocate  a  simpler  procedure  which  appears  to  provide 
the  same  desirable  attributes  to  nearly  the  same  extent.  This  procedure 
consists  of  the  following  two  steps: 

(1)  Use  Prony's  method  on  the  given  data,  but  with  a  prediction  order 
L  which  is  larger  than  the  maximum  number  of  exponentials  which  are  expected 
in  the  signal.  The  result  is  a  set  of  L  exponentials  which  are  candidates 
for  signal  components. 

(2)  Out  of  the  L  exponential  functions  which  are  provided  by  the  high- 
order  Prony  calculation,  determine  the  best  subset  of  size  M.  A  best  subset 
of  M  exponentials  is  one  for  which  a  linear  combination  of  the  M  exponentials 
best  approximates  the  observed  data  using  a  least  squares  criterion.  One 
can  chrck  all  of  the  (£)  possible  subsets  of  size  M  of  the  L  exponential 

n 

functions  to  find  the  best  combination. 

A  simpler  approach  to  step  2  is  to  use  the  procedure  of  Hocking  and 
Leslie  |jxjj  as  we  have  previously  suggested  Dll .  In  the  procedure  of 


6 


Hocking  and  Leslie  a  best  subset  can  usually  be  found  without  searching  over 
all  possible  subsets.  Hocking  and  Leslie  accomplish  this  by  first  solving 
a  related,  but  different,  problem.  They  search  for  the  basis  functions 
(exponentials,  in  our  application)  which  contribute  most  to  the  sunned  magni¬ 
tude-squared  errors  by  their  deletion.  This  provides  an  initial  importance 
ordering  of  the  exponentials.  Hocking  and  Leslie  prove  that  the  fitting 
errors  associated  with  these  single-deletion  sets  provide  convenient  threshold 
values  of  the  error  for  recognizing  the  global  optimality  of  a  particular 
combination  of  M  exponentials  being  tested. 

A 

If  M  is  not  known  apriori,  an  estimate  of  M,  (M)  can  be  found  as  follows. 

A 

Choose  M  *  1,  and  find  the  best  subset  of  size  unity  that  best  fits  the  data. 

A 

Call  the  corresponding  minimum  error,  E^.  Then,  choose  M  ■  2  and  find  the  best 
subset  of  size  two  and  the  corresponding  minimum  error  E^.  Repeat  the  pro¬ 
cedure  until  the  rate  of  decrease  of  the  error  with  increasing  values  of  M 
is  small,  consistent  with  the  modeling  of  broadband  noise.  The  integer  i 
at  which  shows  the  significant  drop  in  rate  of  decrease  is  taken  as  M. 

We  now  give  a  simulation  example. 

III.  Simulation  Results 

If  the  data  is  known  to  be  composed  of  undamped  sinusoids,  as  we  assume 
in  this  example,  forward  and  backward  prediction  equations  can  be  used 
simultaneously  to  obtain  extra  prediction  equations  for  Hildebrand's  least 
squares  form  of  the  Prony  method  [32,33,34j. 

A  sequence  y(n)  consisting  of  two  complex  sinusoids  and  white,  complex 
Gaussian  noise  w(n)  was  generated  using  the  formula  below; 

y<n>  -  .^S”  *  V*  *  *2>  ♦  *(.),  .  -  0,1,2.. .21 

Here  "  *2  "  W1  "  2k(0.52),  ■  2tt ( 0 . 5 )  and  j  ■  -i.  The  variance  of 


v\-»  -I. 


7 


2  2  2 

the  real  or  imaginary  part  of  w(n)  ia  c  .  SNR  is  defined  as  10  log^Q  (a^/2 o  ). 
The  coefficients  of  the  polynomial  B(z)  were  found  by  solving 
the  forward-backward  linear  prediction  equations  as  in  ref.  Q35J  in  the 
least-squares  sense.  L  was  chosen  to  be  12  (N/2).  The  12  zeros  of  G(z) 
were  found  and  the  best  subset  of  2  out  of  the  12  which  minimized  E  in 
formula  11  was  computed.  The  frequency  estimates  of  the  two  sine waves, 

A  A 

and  f2  are  the  angles  of  the  two  chosen  exponents  (divided  by  2tt).  This 
simulation  was  repeated  500  times  and  the  root  mean  square  (RMS)  value  of  the 
frequency  estimation  error  was  computed  at  SNR  values  in  the  range  of  30dB 
to  7dB.  They  are  given  in  Table  1  along  with  the  appropriate  CR  bounds  and 
SVD-method  values  which  were  taken  from  ref.  (35).  Comparing  these  figures 
with  those  in  £.35],  we  note  that  the  SVD-based  methods  are  slightly  better 
in  performance.  This  difference  is  due  to  the  signal  enhancement  achieved  by 


SVD. 


Figure  1  shows  the  minimum  subset  error  Ea  for  different  choices  of  M  at 

”  N  ,  .2 


different  SNR  values.  The  value  E.  at  M  *  0  is  the  data  "energy"  £  y(n) 

A  °  h«l 

Note  the  clear  drop  in  E  at  M  ■  2. 

IV.  Discussion  and  Conclusions 

Ideally,  to  fit  exponentials  to  a  data  sequence  y(n),  one  has  to  mini¬ 


mize  the  error 


N 

£ 

n«l 


y(n)  -  £ 


k-1 


a  skn 
*ke 


with  respect  to  a^'s  end  s^’s 


simultaneously.  This  is  a  difficult  problem  even  if  the  value  of  M  is  known. 

* 

Instead,  we  find  the  exponents  s(k)  separately  as  is  often  done.  However,  we 
have  made  use  of  the  fact  that,  if  the  data  is  composed  of  exponentials  and 
noise,  overestimating  the  degree  L  (>  M)  of  the  polynomial  B(s)  improves  the 
accuracy  of  the  M  signal-zero  locations.  Subsequently,  we  select  the  M  out 
of  the  L  exponentials  that  best  explain  the  data.  The  new  procedure  extends 
the  threshold  of  the  forward-backward  covariance  method 0*2, 33, 3^]  and  is 
only  slightly  inferior  to  SVD  based  methods  (8,35). 


-*^ 


V'v'Ll*-- 


.J*, 


*  * 


Mean  Square  Error 


0.427xl0~3 

(SVD)  0.403  x  10“3(L  ■  12) 

0.311xl0“3 

0.130xl0-2 

0.984xl0~3 

0.238xl0~2 

0.175xl0“2 

-2 

0.373x10 

(SVD)  0.313  x  10_2  (  L  ■  12) 

0.276xl0~2 

0. 417x1 0“2 

0.311xl0~2 

0.601xl0-2 

_ 1 

0.490x1 0“2 

TABLE  1:  Mean  square  error  of  the  frequency  (f^)  estimation 
error  vs.  s.NR.  CRB  stands  for  the  Cramer-Rao  bound  which  is 
the  lower  bound  on  the  standard  deviation  of  the  frequency 
estimation  error  for  an  unbiased  estimator.  The  bias  in  the 
frequency  estimates  was  insignificant  except  at  SNR  »  7dB. 

Below  7dB  the  mean  square  error  is  large  due  to  the  presence 
of  outliers.  For  the  proposed  subset  selection  method  the 
prediction  order  is  twelve  (L  ■  12).  For  comparison,  two  values 
of  the  error  for  the  SVD  method  are  provided  Qs] . 
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